{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sparse covariance estimation for Gaussian variables\n",
    "\n",
    "A derivative work by Judson Wilson, 5/22/2014.<br>\n",
    "Adapted (with significant improvements and fixes) from the CVX example of the same name, by Joelle Skaf, 4/24/2008.\n",
    "\n",
    "Topic References:\n",
    "\n",
    "* Section 7.1.1, Boyd & Vandenberghe \"Convex Optimization\" \n",
    "\n",
    "## Introduction\n",
    "\n",
    "Suppose $y \\in \\mathbf{\\mbox{R}}^n$ is a Gaussian random variable with zero mean and\n",
    "covariance matrix $R = \\mathbf{\\mbox{E}}[yy^T]$, with sparse inverse $S = R^{-1}$\n",
    "($S_{ij} = 0$ means that $y_i$ and $y_j$ are conditionally independent).\n",
    "We want to estimate the covariance matrix $R$ based on $N$ independent\n",
    "samples $y_1,\\dots,y_N$ drawn from the distribution, and using prior knowledge\n",
    "that $S$ is sparse\n",
    "\n",
    "A good heuristic for estimating $R$ is to solve the problem\n",
    "  $$\\begin{array}{ll}\n",
    "    \\mbox{maximize}   & \\log \\det(S) - \\mbox{tr}(SY) \\\\\n",
    "    \\mbox{subject to} & \\sum_{i=1}^n \\sum_{j=1}^n |S_{ij}| \\le \\alpha \\\\\n",
    "                      & S \\succeq 0,\n",
    "    \\end{array}$$\n",
    "where $Y$ is the sample covariance of $y_1,\\dots,y_N$, and $\\alpha$ is a sparsity\n",
    "parameter to be chosen or tuned.\n",
    "\n",
    "## Generate problem data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cvxpy as cp\n",
    "import numpy as np\n",
    "import scipy as scipy\n",
    "\n",
    "# Fix random number generator so we can repeat the experiment.\n",
    "np.random.seed(0)\n",
    "\n",
    "# Dimension of matrix.\n",
    "n = 10\n",
    "\n",
    "# Number of samples, y_i\n",
    "N = 1000\n",
    "\n",
    "# Create sparse, symmetric PSD matrix S\n",
    "A = np.random.randn(n, n)  # Unit normal gaussian distribution.\n",
    "A[scipy.sparse.rand(n, n, 0.85).todense().nonzero()] = 0  # Sparsen the matrix.\n",
    "Strue = A.dot(A.T) + 0.05 * np.eye(n)  # Force strict pos. def.\n",
    "\n",
    "# Create the covariance matrix associated with S.\n",
    "R = np.linalg.inv(Strue)\n",
    "\n",
    "# Create samples y_i from the distribution with covariance R. \n",
    "y_sample = scipy.linalg.sqrtm(R).dot(np.random.randn(n, N))\n",
    "\n",
    "# Calculate the sample covariance matrix.\n",
    "Y = np.cov(y_sample)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve for several $\\alpha$ values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Completed optimization parameterized by alpha = 10, obj value = -16.167608186713004\n",
      "Completed optimization parameterized by alpha = 2, obj value = -22.545759632606043\n",
      "Completed optimization parameterized by alpha = 1, obj value = -26.989407069609157\n"
     ]
    }
   ],
   "source": [
    "# The alpha values for each attempt at generating a sparse inverse cov. matrix.\n",
    "alphas = [10, 2, 1]\n",
    "\n",
    "# Empty list of result matrixes S\n",
    "Ss = []\n",
    "\n",
    "# Solve the optimization problem for each value of alpha.\n",
    "for alpha in alphas:\n",
    "    # Create a variable that is constrained to the positive semidefinite cone.\n",
    "    S = cp.Variable(shape=(n,n), PSD=True)\n",
    "    \n",
    "    # Form the logdet(S) - tr(SY) objective. Note the use of a set\n",
    "    # comprehension to form a set of the diagonal elements of S*Y, and the\n",
    "    # native sum function, which is compatible with cvxpy, to compute the trace.\n",
    "    # TODO: If a cvxpy trace operator becomes available, use it!\n",
    "    obj = cp.Maximize(cp.log_det(S) - sum([(S*Y)[i, i] for i in range(n)]))\n",
    "    \n",
    "    # Set constraint.\n",
    "    constraints = [cp.sum(cp.abs(S)) <= alpha]\n",
    "    \n",
    "    # Form and solve optimization problem\n",
    "    prob = cp.Problem(obj, constraints)\n",
    "    prob.solve(solver=cp.CVXOPT)\n",
    "    if prob.status != cp.OPTIMAL:\n",
    "        raise Exception('CVXPY Error')\n",
    "\n",
    "    # If the covariance matrix R is desired, here is how it to create it.\n",
    "    R_hat = np.linalg.inv(S.value)\n",
    "    \n",
    "    # Threshold S element values to enforce exact zeros:\n",
    "    S = S.value\n",
    "    S[abs(S) <= 1e-4] = 0\n",
    "\n",
    "    # Store this S in the list of results for later plotting.\n",
    "    Ss += [S]\n",
    "\n",
    "    print('Completed optimization parameterized by alpha = {}, obj value = {}'.format(alpha, obj.value))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Result plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArMAAALACAYAAACAUfp0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3dtxG9mWJuC1JmQAR6f53BGQB5TKgmZ5wFNtwbA8UMWxoILlAXksOC15IHlQpDwQO+aZPSp6sOcBCSkFkcSFABIL+X0RCAlEMnPvvCz+mchLttYCAAAq+l9DNwAAANYlzAIAUJYwCwBAWcIsAABlCbMAAJQlzAIAUJYwCwBAWcIsAABlHXSYzcyTzLzIzJaZN5n5dug27UpmnnV9f5eZk6HbU0FmXmbm5dDt2GfmEduSmZP5et17XWbmX5l59ozx7/26m5lHXc1e+DSjCv05NBXnebddvc3Mi0c+f9vlhfPMPN91+zblxdAN2KbW2qeI+NQVwH+11v4Yuk27kJlHEXHRWnu1qPhn5nlr7WpHTdt374ZuwBBWXAdGOY/YvtbabUT89li97kLE6TLjemSd3sm6+5ya2lq7j4i/LxNmw7a4EYdc/zLzNCKOIuLVI59fRMSfrbX3s/eZeTZ7X8lBh9kRexMRtxERS6yUr7ffnBpaax+HbsNAll4HRjyPGFhr7VNmvlly8B/W6R2uuzupqbbFjTnY+jdrb2b+FNNQO++8tfZb7/2HiPgtIsqF2YM+zYCndaddvBy6HQzHOsC+m/t26XqJ4Qdbp21PtYx5eWXmyQM//hJLfvuxb0Z5ZLY79H4R06OXs/Nffo6Iz621q654zs4v+Xt3ROAsIv4ZEVettd+6r/L/ERF/RsRPEfGhtfaxN+7rmO7l/GdMvzJ73204n2K6h/RzTE8FuH1sXE+0/ygizrv2R0RMZl/JdSvo3yNi0k3vY3e6xUPz4KeIOOmGu+/6/mD7I+K++/mX1trPXRv+GRFnEfG/u6/HYs2+zIa/j4ij3lceD/Zz0fLp2js7T/i7NjzRv9vu97+01n7uzctF43lwHVqhf0vNq0emd9KN82P3/5cR8bq19mvv9xb1YZV14Lt51M3zf3Tj/49uOdzEdP2+cPoKz9VtI/8Z3ZGifi17qJ7GdF18aJ0+ie/X3Z1tT71+PLqtd1/3fo5pmLhfYr4s05+V/qYtmuYDbdh47X6sHWOof3PzbBIRn7Z8JPhlTNe3voXr3t5qrR38K6ZF4u3cz866n0+690fT2fH189OIuJn7nfO5cR713t/M3sd0hfzcjfOke51HxNnc9E8WjeuR/tzMDX8WEZdzbf+wxHw5i4h3D/z8h/b3hv8wN2yba8s6fTmZ+/2Th36338+nlk/3e/15/deS/ftuvi0xnifXoSX6t+q8enB6EXHaG+bDXJuX6cMq68AP69Zcny7W3U69vFr7uj59iIi3MQ0Vfz0wzFP19LF1en773uX29NTfiw9zNWIyX0cemU9L9Wdu+Ef/pq2xnDZeuxdM72DrX7fM5+fZbPqTXp8vl3wdzY3/InoZoT8/5342m6eP/h3a19coj8x27iO+XnQQrbX7zPz6YZseZZ1k5kmb7nGdRsR/RXz92uu+dUcjO9cR8UtMjwx+iYjb7vNP3e+8jIjL7t+P7dse7KJxfWf2lVt/+DY96vsuM3+bG8+6fmh/58lxr9GX05hupP1p9PfWn+rno8unG8dtb5yRmUe9cT3Wv3mLxvPkOrRE/5aeVwumdzs3TP9rs0V9eMyy8yhiegToQ3d06fcFw8IyPrRv3zb9vwc+v40H6umKdrI9PbWtZ+Z1zNWINv22bo3urP83bVVbrN3r9O8Q6t+7iPh9rl33Mb3+ZdK6o7ytd9R5A+bnVTzwvowxh9mI7zeCiPhhRf89In7tXpP27ZD/pBu2f27Ju/hxo/qq24Bn47rMzE8R8R9LjqtvEj9+NTCb3iQWb3jLWicUr9OX7z7rFcdl+vnY8vky+9qoN46X8X2flunfMuN5ah1a1L9V5tWj04vv59P8PFumD49Zah3o/vi+j+kfDqcWsGlfg+rsSuvH6ukaO/O72J6e2tZPF/zuqtb9m7aqbdXuRQ6u/uX01pknD+yQ3cY0KP9r1XEu6Uv8eFHYUcT3OyJVjD3MLnIVEf/d7XH1N6LbmB6GX7oQZOZpN/zsfJ2L+HZ+zCrjuo2H956OYnEQeqp9k4hve71rWqcvj90Dd5l+PrZ8bqI7ShARkZn/7P5dZo+877njWdS/ldahNS3dh3XXge5crz+jO0+7jeQWeOzG3Pr4NRg+Uk/nb+e1ibrWt/L2FE9s65n5VI3Yhsdq5qq2Vbs3rUL9O4knDkK1788TX/Yetwu/pe2Ois8P8zK6baqasd/N4MlD6t3KcB3Tc00+9n7+PiLue8VqdmPih64OnDmZ2zP/1zrj6oY/mhv+LCLer3lUYjaeyRIb8XcFbL6Na/TlY0Tc9udL7+unhf18aPl003rZK179Pc+lr9JcYTyPrkNL9G/VdejJ6c1/vmQfVl0HvtON8x9df/5PRPy6RB9gZd269lP39sF62lllnd769vTUtt6rESe9z56z/az1N63XpqUeSrGN2r2kUvWvmyeL/u7MLmLsj3fStetz/+ettV+XfC2bBf5rbpn/HN8usKul7cGJu9t6xXSP5yKmJzTfRHcRWPfzDxHxV0wvMjjqDXcZ3cnW3bCnMXfxWPt2ovRFTE+iPovuJPRu+O/G3f38vHvNhn+7aFxP9Ks//PncuPp9u4gFJ3LH9Ouu8/h28dSD7e8NPzsCctZN669uHJNn9OWyN87TZfr51PLpxve2++y0+/2Lrr2PLZ/+fDtfYjxLrUMr9O/JebXM9OLbRQs38W19fLQPq64D8/Oo++xzfH8RyufZ7w29/XvVenXr8Hf1uve67Nar2UVEj9bT7vP5dXp+3d3Z9tT97NFtfe6z2XhbN54H6/c6/en97mN/097GEhcOz7V7o7X7iWmVrH/dMD9cwPhA/2br86yNs7swXPbXozW2qVkG+ty9vut7r41n3b9rT2voV3adAQBgg3qnxLBFYz/NAABgW+YvsmILhFkAgA3rzn3d5gVudJxmAABAWY7MAgBQ1l7cZ3b2hJSYXmV3sPeo7G7jMbtNx09tjedhV5SZF4fe1+w9A7yt9zSiUnrb7Nen0zAeavZhG0PNjhhX3T70mj34kdneo+8+xvRefEvfC7SgX2J637vZo2zPB27P1nXLc5c3BR/Kr91ynfTvr3iIumV62x64RyaHT80+bCOq2REjqdtjqNmDh9mY3gR7doL0bUzvi3aQWmtXvT2iSRR90sayxnLye/cH7iYzJ621P9rmnja0r64j4t3sqEb7/vnsHD41+0CNpWZHjK5uH3zN3ocwO3/bir8N0ood6grGlwPfeCLWeKJKUa+615fMvJx70szBadOny1zG9GbjrwZuDrunZh+usdTsiBHV7THU7H0Is/ex+BF1h+astfbr0I3YphHeKPpzVzBuYvokl4PVfc38sbX2qvee8VCzD9AIa3bESOr2GGr2PoTZP+Pbnv4kpo+LO1iZeTa7YOLAzzX7kpmn3UYzOcRzdHr+7P3/KKZ/7A9Z/2uq32N8wWbs1OzDNKaaHTGuun3wNXvwMNs7+fq0e3+we4ZdHy8y8yYzb4Zuzza11j51y/JlHPgTULp1+Ki3Dh/claJzrjLzvOvvLyPoLz1q9mEaU82OGF3dPvia7aEJAACUNfiRWQAAWJcwCwBAWcIsAABlCbMAAJQlzAIAUNbehNkxPPN6Zkx9jRhXf/WVsRjb8h9Tf/X1MB1yX/cmzMYBP33jAWPqa8S4+quvjMXYlv+Y+quvh+lg+7pPYRYAAFay0YcmZOZonsDw+vXrtX/37u4ujo+PN9ia/Tam/urrcm5uBnuY0v+01saxgJagZi/Ptn2YKvZ1iPr53O1nXTc3N0vVbGF2TZ6cBuvLzKEmfdNaezPUxPeNmg31DFE/h9p+MnOpmu00AwAAyhJmAQAoS5gFAKAsYRYAgLKEWQAAyhJmAQAoS5gFAKAsYRYAgLKEWQAAynqxaIDMPIuI+4g4aa39sf0mAbAuNRsYmyePzHZFMVprHyPiPjNPd9IqAFamZgNjtOg0g58i4rb7/21EnGy3OQA8g5oNjM6iMHs09/5v22oIAM+mZgOjs+ic2fuIePnUAJl5HhHnG2sRAOtSs4HRWRRm/4xve/qTiPgwP0Br7SoiriIiMrNttHUArELNBkbnydMMWmvvI2Iyu4igu6gAgD2kZgNjlK1tbsd8THv5m5xvMDaZOdSkb1prb4aa+L5Rs6GeIernUNtPZi5Vsz00AQCAsoRZAADKEmYBAChLmAUAoCxhFgCAsoRZAADKEmYBAChLmAUAoCxhFgCAsoRZAADKejF0A6oa8HGcg/AoyMM1tnUZdmmo7UvNZkwcmQUAoCxhFgCAsoRZAADKEmYBAChLmAUAoCxhFgCAsoRZAADKEmYBAChLmAUAoCxhFgCAsoRZAADKEmYBAChLmAUAoCxhFgCAsoRZAADKEmYBAChLmAUAoCxhFgCAsoRZAADKEmYBAChLmAUAoCxhFgCAsoRZAADKEmYBAChLmAUAoCxhFgCAsoRZAADKEmYBACjrxSZH9vr167i+vt7kKJeSmTuf5lBaa0M3gQMzpnVqTLViGUPV7CEMtezHtH3BUByZBQCgLGEWAICyhFkAAMoSZgEAKEuYBQCgLGEWAICyhFkAAMoSZgEAKEuYBQCgLGEWAICyhFkAAMoSZgEAKOvFUx9m5lFEnHZvf2qt/bb9JgGwDjUbGKNFR2Z/iYiXrbX3ERGZeb79JgGwJjUbGJ0nj8y21q56bycRcbnd5gCwLjUbGKOlzpnNzElEfGmt3T7w2XlmXmfm9d3d3cYbCMBq1GxgTJa9AOystfbrQx+01q5aa29aa2+Oj4832DQA1qRmA6OxMMxm5llr7Y/u/6eLhgdgOGo2MDZPhtmuEF5k5k1m3uyoTQCsQc0GxmjRBWAfI+LVjtoCwDOo2cAYeWgCAABlCbMAAJQlzAIAUJYwCwBAWcIsAABlCbMAAJQlzAIAUJYwCwBAWcIsAABlPfkEMB7XWhtkupk5yHSH6i/bN8Q6ZX0atzGtc2r24Rpq2fIjR2YBAChLmAUAoCxhFgCAsoRZAADKEmYBAChLmAUAoCxhFgCAsoRZAADKEmYBAChLmAUAoCxhFgCAsoRZAADKEmYBAChLmAUAoCxhFgCAsoRZAADKEmYBAChLmAUAoCxhFgCAsoRZAADKEmYBAChLmAUAoCxhFgCAsoRZAADKEmYBAChLmAUAoCxhFgCAsoRZAADKEmYBACjrxdAN2ITW2tBN2Jmh+pqZO5/mmPoaMVx/x7T9sB/GtM6NqY6Nqa8R41qP950jswAAlCXMAgBQljALAEBZwiwAAGUJswAAlCXMAgBQljALAEBZwiwAAGUJswAAlCXMAgBQljALAEBZS4fZzLzYZkMA2Bw1GxiLpcJsZp5GxGTLbQFgA9RsYEwWhtnMnETE7Q7aAsAzqdnA2CxzZHbSWnu0MGbmeWZeZ+b13d3dBpsGwBrUbGBUngyzmXnaWvv41DCttavW2pvW2pvj4+PNtg6ApanZwBi9WPD5l+7cq6OImGTmSWvt0w7aBcDq1GxgdJ48Mtta+9Tt5b+MaXEEYE+p2cAYLXU3g+5rqVf28AH2n5oNjImHJgAAUJYwCwBAWcIsAABlCbMAAJQlzAIAUJYwCwBAWcIsAABlCbMAAJQlzAIAUNaLoRtADa21nU8zM3c+zYhh+gqwSWo2Y+LILAAAZQmzAACUJcwCAFCWMAsAQFnCLAAAZQmzAACUJcwCAFCWMAsAQFnCLAAAZQmzAACUJcwCAFCWMAsAQFnCLAAAZQmzAACUJcwCAFCWMAsAQFnCLAAAZQmzAACUJcwCAFCWMAsAQFnCLAAAZQmzAACUJcwCAFCWMAsAQFnCLAAAZQmzAACUJcwCAFCWMAsAQFnCLAAAZb0YugHUkJk7n2ZrbefTjBimrxHD9Rc4PGr29qnZ+8ORWQAAyhJmAQAoS5gFAKAsYRYAgLKEWQAAyhJmAQAoS5gFAKAsYRYAgLKEWQAAyhJmAQAoS5gFAKCsF4sGyMyTiJhERLTW3m+9RQCsTc0GxmaZI7O/dgVxkpmTbTcIgGdRs4FRefLIbGaeR8RNZk5aa3/sqE0ArEHNBsZo0ZHZV93rS2ZeZubRDtoEwHrUbGB0ljnN4HNr7T4ibiLifP7DzDzPzOvMvL67u9t4AwFYiZoNjMqiMPtn7/9HEXE/P0Br7aq19qa19ub4+HijjQNgJWo2MDpPhtnuIoKjzDzt3l/tpFUArEzNBsZo4a25ehcRfNxyWwB4JjUbGBsPTQAAoCxhFgCAsoRZAADKEmYBAChLmAUAoCxhFgCAsoRZAADKEmYBAChLmAUAoCxhFgCAshY+zhYiIlprQzdhZ4bqa2YOMt0xLVvYNdv19qnZODILAEBZwiwAAGUJswAAlCXMAgBQljALAEBZwiwAAGUJswAAlCXMAgBQljALAEBZwiwAAGUJswAAlCXMAgBQljALAEBZwiwAAGUJswAAlCXMAgBQljALAEBZwiwAAGUJswAAlCXMAgBQljALAEBZwiwAAGUJswAAlCXMAgBQljALAEBZwiwAAGUJswAAlCXMAgBQljALAEBZL4ZuwCZk5s6n2Vrb+TQjhulrxHD9HZMxrVPWp/1wc3MzWE2B6tTs/eHILAAAZQmzAACUJcwCAFCWMAsAQFnCLAAAZQmzAACUJcwCAFCWMAsAQFnCLAAAZQmzAACUJcwCAFDWi0UDZOZZRNxHxKS1drX9JgGwLjUbGJsnj8xm5mlE3LbWPkbEbWae7KZZAKxKzQbGaNFpBtcR8a4riJPW2qcdtAmA9ajZwOg8GWZba/cRcRkR7yLi1UPDZOZ5Zl5n5vXd3d0WmgjAMlat2TttHMCWLDrN4CwiPrbWXvXef6e1dtVae9Nae3N8fLylZgKwyKo1e+cNBNiCRacZ9L+m+j0iXm65PQCsT80GRmfR3QyuMvM8Im7DlbEA+07NBkbnyTDbnX+lGAIUoGYDY+ShCQAAlCXMAgBQljALAEBZwiwAAGUJswAAlCXMAgBQljALAEBZwiwAAGUJswAAlLXocbYrubm5iczc5CiZ01obugkcmCHWKXViP7x+/Tqur6+HbgawAjX7R47MAgBQljALAEBZwiwAAGUJswAAlCXMAgBQljALAEBZwiwAAGUJswAAlCXMAgBQljALAEBZwiwAAGUJswAAlCXMAgBQljALAEBZwiwAAGUJswAAlCXMAgBQljALAEBZwiwAAGUJswAAlCXMAgBQljALAEBZwiwAAGUJswAAlCXMAgBQljALAEBZwiwAAGUJswAAlCXMAgBQVrbWNjeyzLuI+L9r/vq/RcT/bKwx+21MfY0YV3/1db/9e2vteOhG7As1eyVj6q++HqaKfV2qZm80zD5HZl631t4M3Y5dGFNfI8bVX31lLMa2/MfUX309TIfcV6cZAABQljALAEBZ+xRmr4ZuwA6Nqa8R4+qvvjIWY1v+Y+qvvh6mg+3r3pwzCwAAq9qnI7MAALASYRYAgLKEWQAAyhJmAQAoS5gFAKAsYRYAgLKEWQAAyhJmAQAoS5gFAKAsYRYAgLKEWQAAyhJmAQAoS5gFAKAsYRYAgLKEWQAAyhJmAQAoS5gFAKAsYRYAgLKEWQAAyhJmAQAoS5gFAKAsYRYAgLJeDN2AXcnMSUT8GhFvI+JTRPyr9/GriPglIv5Pa+39muO/jIhorf36zKZuTWYeRcQ/I+KstZYLht37/hyaavO8W5/Ou7c/RcSH1trVgE3igKjZava+qzjPu+3qLCL+1lr7bej2bMpowmxr7TYifsvMs4j4V2vtj/7n3Up5usy4MvP8gT/a7zbT0rWmvZTW2n1E/D0z2xKD76Q/h27F5VVtnv+jXwwz83NmhkDLJqjZavYQDrlmZ+ZpRBzFdGfwoIwmzC7SWvuUmW+WHPz1A7//ccNNWnra27DD/hy6pZdXpXneHTGazP34MiJ+iwhhlq1Ts79XqX7suYOs2RHf2puZP8U01B6M0Z8z2+31z1wvMfzbiHi5vRbt57RZ3YEvr5cRcdp9ZTVzHz8GXNgoNZttsbzqGvWR2e7o0n9GxPuI6Z5+77PZeVpHEfFzRFzE9A/1TxFx0n1+31q7ysyTmJ7X9KW19nN3KP8iIm5jerQqIuIkpn/sP3b/fxkRr/vn2nTjmYWB2TmIsz2p04em3evHPyLiz/nf6z6/iIjPEfGla8Oi+bJMf36OiM9d/8+6zyMi/t4dMTnrxnG1znk5c326j4ij2blxvXM1b7vBJ621P9Ztxy6XV29a1xHxIabr37+6affn+VnX/0lE/EfXl5uYro8Xq35tOTfPJhHx6TlHFbqvgP/33I9/jun8gq1Qsx+dL2q2mj1urbVRvWJaID7E9KKCi4j464FhzmN6wv3s/VlEnPT+/+6B3zmN6YbR/53PMd1oI6YrdIuI094wH+amczP3/q+5aTw27c8xLRz98Rz1pnHS+2wyXewL59NS/Zkb/mZ+Pj5jOd3Mtftzbxl87V+vbZfPaceOl9f5bJnFtOie9Nr+YW7Yfr8v1pyXkwfm2Wz6k16fL5d8HT0wjaOI+Ku/zLy8NvEKNVvNfnh6avaaNTum29Hlust6H19jPTL7oXUXE2Tm/3vg89uIuMzMlxHxsa13tex9xNcjWNFau8/M2bj7w/S/0vj7bPiZzDxq04sAHtTtDd7PDXMdEb9k5nVMV/qvRy9aa7ddOzbVn9l4P2bmJDNP2nSP9DQi/mudCXW/+1274/u99uj3t7X2PjPfZeZvz2jHTpZX50tE3HbDfVow7M8R8aE7UvP7wl487F1E/D7XrvuIeBPT+XwV8ewrcv8Z3dGIZ4wDHqNmb64/s/Gq2THqmn1QRn/ObHRfV0V8OxerTQ/l/xrd1zKZedMd8l/V7QM/+/LI/yMivmTm28w8750Xtuj8nUnEtJjMXjHdED7GdMVf+BXVCn7oz9x8+T2m8y1iusGtO+3J/LR6RXISP863iO/P11y3HbtYXjNLtakrvO9jWoRXnp/dOa0nD/xxv43p+r3wnMMlpvE2pnv5giy7oGYvT82e+0zNPkxjPTL71dxe2tci0xXH2bk0FzH9mmH+1jCTB8bxHDfR7c124/9n9+8Pe4757cKb25h+hfDDuTSZOTvXZleuIuK/u/n1nHnyVLtv4+Hic9Sb5qbascjKy2vVdaX7w/NnREwy822buz3REk7iiSMJ7ftzDi8fG27Ob7P+dX8Qvp7L1dt2YCvU7I1Ss2NcNftQjT7MznQr4E/d25Oc3i9zVmz+Fd/uZ9jfaCcL/nAv2uP7+nl3YvrL3kbW33s+jeme3kPTvs3Mf2TmZLbRdRvhUff1ze3s65vedNb1ZH+6r3muY3qU7uf+Z0/sbT40nlm7T3shadan95l5Mdffs4h4P9tYn2rHc/oXm1leS+vG+Y/W2m+Z+TEibjLzY2+as3Onnhrv7IKY/ngnXbs+9H++6ldW3RGllxHxsfeH2kVg7ISavRQ1+1t/Rl+zD1rbgxN3d/GK6YpwEdMTxG9iejHB7HUZ04tXZiekn3evs+71dm5c72bDdO9PYrqS/dX9vP/+bXRXM3bTvuzaMjuh/Ca6E9a7z97GdMM67aZ9Ed+fVP/dtLufzcY/a+/pI5/Nxtu68fxwIc+6/en97un8/Op+/jbmTpRfsLyOunHPlsNjfTp/ZHoPtmOJ/m51eXXDfTetR+b52256/Qs6Ps9+rzdPf7gY5oH+zdbnWRtnV/RexpoXfMS3iy3mXz9cPOHltc4r1Gw1+/FpqdnrbVMnXX8/d6/v+l75lV0HgYJ8rQ9Qh5q9HS4Ag9oO6ikuAAdOzd4CYRaK6s6j2ubFEgBsiJq9PU4zAACgLEdmAQAoay9uzdXdpuM+plfVrXpPtjK623bMbhfzU1vj+dcVZebFofc1e8/8bus9faiU3jb79Wk0jIeafdjGULMjxlW3D71mD35kdu4JLvfdvSsP1S8xvc/d+4iIzDwfuD1b1y3PXd4EfCi/dst10rvn6kHqlult+3bPzOfcB5Ni1OzDNqKaHTGSuj2Gmj14mI3pTa9nJ0TfxvQ+aAeptXbV2yOaxIHfXH4sJ7t3f+BuupuC/9E293ShfXUdEe9mRzWax9iOjZp9oMZSsyNGV7cPvmbvQ5idv03F3wZpxQ51BePLgW88EdON5tD7GBHxqnt9yczLXO+Z8GW06VN7LmN6c/FXAzeH3VOzD9dYanbEiOr2GGr2PoTZ+1j8SLpDc9YO/DF0I7wx9OeuYNzE9MktB6v7mvlja+1V7z3joWYfoBHW7IiR1O0x1Ox9CLN/xrc9/R+ePXxoMvNsdsHEgZ9r9iUzT7uNZnKI5+j0/Nn7/1FM/9gfsv7XVL/H+ILN2KnZh2lMNTtiXHX74Gv24GG2d/L1aff+YPcMuz5eZOZNZt4M3Z5taq196pblyzjwJ5506/BRbx0+uCtF51xl5nnX319G0F961OzDNKaaHTG6un3wNdtDEwAAKGvwI7MAALAuYRYAgLKEWQAAyhJmAQAoS5gFAKCsvQmzY3jm9cyY+hoxrv7qK2MxtuU/pv7q62E65L7uTZiNA376xgPG1NeIcfVXXxmLsS3/MfVXXw/TwfZ1n8IsAACsZKMPTcjMQZ7A8Pr16yEmu7a7u7s4Pj4euhk7M6b+6ut+u7m5+Z/WWq1Gb5GavbyK6/u69PUwVezrsjX7IMKsp5gBy8jMm9bam6HbsS/UbGCfLVuznWYAAEBZwiwAAGUJswAAlCXMAgBQljALAEBZwiwAAGUJswAAlCXMAgBQljALAEBZLxYNkJlnEXEfESettT+23yQA1qVmA2Pz5JHZrihGa+1jRNxn5ulOWgXAytRsYIwWnWbwU0Tcdv+/jYiT7TYHgGdQs4HRWRRmj+be/21bDQHg2dRsYHQWnTN7HxEvnxogM88j4nxjLQJgXWo2MDqLwuyf8W1PfxIRH+YHaK2QuRFtAAAJpklEQVRdRcRVRERmto22DoBVqNnA6Dx5mkFr7X1ETGYXEXQXFQCwh9RsYIyytc3tmA+1l7/JPgCHKzNvWmtvhm7HvlCzgX22bM320AQAAMoSZgEAKEuYBQCgLGEWAICyhFkAAMoSZgEAKEuYBQCgLGEWAICyhFkAAMoSZgEAKOvFJkf2+vXruL6+3uQol5KZO5+mxzEC1anZwCFwZBYAgLKEWQAAyhJmAQAoS5gFAKAsYRYAgLKEWQAAyhJmAQAoS5gFAKAsYRYAgLKEWQAAyhJmAQAoS5gFAKAsYRYAgLKEWQAAyhJmAQAoS5gFAKAsYRYAgLKEWQAAyhJmAQAoS5gFAKAsYRYAgLKEWQAAyhJmAQAoS5gFAKAsYRYAgLKEWQAAyhJmAQAoS5gFAKCsF0M3YBNaazufZmbufJoRw/QVYJPUbGCTHJkFAKAsYRYAgLKEWQAAyhJmAQAoS5gFAKAsYRYAgLKEWQAAyhJmAQAoS5gFAKAsYRYAgLKEWQAAynrx1IeZeRQRp93bn1prv22/SQCsQ80GxmjRkdlfIuJla+19RERmnm+/SQCsSc0GRufJI7Ottave20lEXG63OQCsS80Gxmipc2YzcxIRX1prt1tuDwDPpGYDY7LsBWBnrbVfH/ogM88z8zozr+/u7jbYNADWpGYDo7EwzGbmWWvtj+7/p/Oft9auWmtvWmtvjo+Pt9FGAJakZgNj82SY7QrhRWbeZObNjtoEwBrUbGCMFl0A9jEiXu2oLQA8g5oNjJGHJgAAUJYwCwBAWcIsAABlCbMAAJQlzAIAUJYwCwBAWcIsAABlCbMAAJQlzAIAUJYwCwBAWU8+zpbHtdYGmW5mDjLdofrL9g2xTlmf2DU1Gw6XI7MAAJQlzAIAUJYwCwBAWcIsAABlCbMAAJQlzAIAUJYwCwBAWcIsAABlCbMAAJQlzAIAUJYwCwBAWcIsAABlCbMAAJQlzAIAUJYwCwBAWcIsAABlCbMAAJQlzAIAUJYwCwBAWcIsAABlCbMAAJQlzAIAUJYwCwBAWcIsAABlCbMAAJQlzAIAUJYwCwBAWcIsAABlCbMAAJT1YugGsJrW2iDTzcydT3Oovo6N+Qzbo2bD9jkyCwBAWcIsAABlCbMAAJQlzAIAUJYwCwBAWcIsAABlCbMAAJQlzAIAUJYwCwBAWcIsAABlCbMAAJS1dJjNzIttNgSAzVGzgbFYKsxm5mlETLbcFgA2QM0GxmRhmM3MSUTc7qAtADyTmg2MzTJHZiettUcLY2aeZ+Z1Zl7f3d1tsGkArEHNBkblyTCbmaettY9PDdNau2qtvWmtvTk+Pt5s6wBYmpoNjNGLBZ9/6c69OoqISWaetNY+7aBdAKxOzQZG58kjs621T91e/suYFkcA9pSaDYzRUncz6L6WemUPH2D/qdnAmHhoAgAAZQmzAACUJcwCAFCWMAsAQFnCLAAAZQmzAACUJcwCAFCWMAsAQFnCLAAAZb0YugHU0Frb+TQzc+fTjBimrwCbpGYzJo7MAgBQljALAEBZwiwAAGUJswAAlCXMAgBQljALAEBZwiwAAGUJswAAlCXMAgBQljALAEBZwiwAAGUJswAAlCXMAgBQljALAEBZwiwAAGUJswAAlCXMAgBQljALAEBZwiwAAGUJswAAlCXMAgBQljALAEBZwiwAAGUJswAAlCXMAgBQljALAEBZwiwAAGUJswAAlCXMAgBQ1ouhGwCPaa0NMt3MHGS6Q/UXYBPUbIbiyCwAAGUJswAAlCXMAgBQljALAEBZwiwAAGUJswAAlCXMAgBQljALAEBZwiwAAGUJswAAlCXMAgBQ1otFA2TmSURMIiJaa++33iIA1qZmA2OzzJHZX7uCOMnMybYbBMCzqNnAqDx5ZDYzzyPiJjMnrbU/dtQmANagZgNjtOjI7Kvu9SUzLzPzaAdtAmA9ajYwOsucZvC5tXYfETcRcT7/YWaeZ+Z1Zl7f3d1tvIEArETNBkZlUZj9s/f/o4i4nx+gtXbVWnvTWntzfHy80cYBsBI1GxidJ8NsdxHBUWaedu+vdtIqAFamZgNjtPDWXL2LCD5uuS0APJOaDYyNhyYAAFCWMAsAQFnCLAAAZQmzAACUJcwCAFCWMAsAQFnCLAAAZQmzAACUJcwCAFCWMAsAQFkLH2cLY9NaG2S6mTnIdIfqL8AmqNk4MgsAQFnCLAAAZQmzAACUJcwCAFCWMAsAQFnCLAAAZQmzAACUJcwCAFCWMAsAQFnCLAAAZQmzAACUJcwCAFCWMAsAQFnCLAAAZQmzAACUJcwCAFCWMAsAQFnCLAAAZQmzAACUJcwCAFCWMAsAQFnCLAAAZQmzAACUJcwCAFCWMAsAQFnCLAAAZQmzAACUJcwCAFCWMAsAQFkvhm4Aq8nMQabbWhtkumMy1DweYp2yPgHVqdn7w5FZAADKEmYBAChLmAUAoCxhFgCAsoRZAADKEmYBAChLmAUAoCxhFgCAsoRZAADKEmYBAChr4eNsM/MsIu4jYtJau9p+kwBYl5oNjM2TR2Yz8zQibltrHyPiNjNPdtMsAFalZgNjtOg0g+uIeNcVxElr7dMO2gTAetRsYHSeDLOttfuIuIyIdxHxaictAmAtajYwRotOMziLiI+ttVe99/PDnGfmdWZe393dbamZACyiZgNjtOg0g/7XVL9HxMv5AVprV621N621N8fHxxtvIABLU7OB0Vl0N4OrzDyPiNtwZSzAvlOzgdF5Msx2518phgAFqNnAGHloAgAAZQmzAACUJcwCAFCWMAsAQFnCLAAAZQmzAACUJcwCAFCWMAsAQFnCLAAAZQmzAACU9eTjbNk/rbWhm8CBGWKdysydTxPgEKjZP3JkFgCAsoRZAADKEmYBAChLmAUAoCxhFgCAsoRZAADKEmYBAChLmAUAoCxhFgCAsoRZAADKEmYBAChLmAUAoCxhFgCAsoRZAADKEmYBAChLmAUAoCxhFgCAsoRZAADKEmYBAChLmAUAoCxhFgCAsoRZAADKEmYBAChLmAUAoCxhFgCAsoRZAADKEmYBAChLmAUAoCxhFgCAsrK1trmRZd5FxP9d89f/LSL+Z2ON2W9j6mvEuPqrr/vt31trx0M3Yl+o2SsZU3/19TBV7OtSNXujYfY5MvO6tfZm6Hbswpj6GjGu/uorYzG25T+m/urrYTrkvjrNAACAsoRZAADK2qcwezV0A3ZoTH2NGFd/9ZWxGNvyH1N/9fUwHWxf9+acWQAAWNU+HZkFAICVCLMAAJQlzAIAUJYwCwBAWcIsAABl/X/XGVabMuWSXwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x864 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Show plot inline in ipython.\n",
    "%matplotlib inline\n",
    "\n",
    "# Plot properties.\n",
    "plt.rc('text', usetex=True)\n",
    "plt.rc('font', family='serif')\n",
    "\n",
    "# Create figure.\n",
    "plt.figure()\n",
    "plt.figure(figsize=(12, 12))\n",
    "\n",
    "# Plot sparsity pattern for the true covariance matrix.\n",
    "plt.subplot(2, 2, 1)\n",
    "plt.spy(Strue)\n",
    "plt.title('Inverse of true covariance matrix', fontsize=16)\n",
    "\n",
    "# Plot sparsity pattern for each result, corresponding to a specific alpha.\n",
    "for i in range(len(alphas)):\n",
    "    plt.subplot(2, 2, 2+i)\n",
    "    plt.spy(Ss[i])\n",
    "    plt.title('Estimated inv. cov matrix, $\\\\alpha$={}'.format(alphas[i]), fontsize=16)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
